function [Ytest d C] = mdm_freq(COVtest,COVtrain,Ytrain,varargin)

    if isempty(varargin)
        method_mean = 'riemann';
        method_dist = 'riemann';
    else
        method_mean = varargin{1};
        method_dist = varargin{2};
    end
    
    labels = unique(Ytrain);
    Nclass = length(labels);
    Nfreq =  size(COVtrain,4);
    C = zeros(size(COVtrain,1),size(COVtrain,1),Nclass,Nfreq);
    NTesttrial = size(COVtest,3);
    
    %% training - estimation of center
    for i=1:Nclass
        
        %get all trials for a specific class
        X = COVtrain(:,:,Ytrain==labels(i),:);
        trialsForThisClass = size(X,3);
        
        % g is all cov trials for specific class, specific frequency
        g = zeros(size(COVtrain,1),size(COVtrain,1),trialsForThisClass);
        
        for f=1:Nfreq
            for t=1:trialsForThisClass
            
                  chunk = X(:,:,t,f);
                  covm = covariances(chunk);

                  g(:,:,t) = covm; %
%                   if (t==1)
%                       C(:,:,i,f) = covm;
%                   else
%                       g = zeros(size(covm,1),size(covm,2),2);
%                       g(:,:,1) = C(:,:,i,f); % current mean
%                       g(:,:,2) = covm; % current input cov (that is to be added to the mean)
%                       %C(:,:,i,f) = mean_covariances(g,method_mean);
%                       C(:,:,i,f) = geodesic(g(:,:,1),g(:,:,2),1/t,method_mean);
%                   end;  
            end;
            
            C(:,:,i,f) = mean_covariances(g,method_mean);
        end;       
    end

    %% classification
    
    % contains the distances for each class i and the frequencies
    % each row is trial
    % each column is the combined distnace to all frequencies for a
    % specific class
    d = zeros(NTesttrial,Nclass);
    
    for j=1:NTesttrial
        for i=1:Nclass
          for f=1:Nfreq
              
            test = covariances(COVtest(:,:,j,f));
            d(j,i) = d(j,i) + distance(test,C(:,:,i,f),method_dist);
            
          end;
        end
    end
    
    [~,ix] = min(d,[],2);
    Ytest = labels(ix);